#version 120
#extension GL_EXT_gpu_shader4 : enable //access to integer and bit operation
uniform mat4 invViewMatrix;
uniform vec3 translateVec3;

uniform usampler3D nodePool3D;
uniform sampler3D blockPool;

uniform sampler1D tranf;

uniform float aspectRatio;

uniform float stepsize;
uniform float isoValue;
uniform float var;
uniform float integralMethod;


//==view
uniform float insideView;

//=====Parameter
const float voxelMax = 255.0;

uniform uint blockPoolDim;
uniform uint nodePoolDim;
uniform uint blockSize;
//uniform float 
float invBlockPoolDim = 1.0/float(blockPoolDim);
///////////////////////////////////////////////////////////////////////
//record rays

//===================================
// structures
struct nodeTexel
{
	uint first;
	uint second;
};

struct ray
{
	vec3 o;
	vec3 d;
};

struct range
{
	float tNear;
	float tFar;
};

struct BlockBox
{
	vec3 boxMin;
	vec3 boxMax;
	bool hasBlock;
	uvec3 blockPointer;
	float sampleRate;
	float LODSamplingFactor;
	float singleColor;
};

struct sample
{
	vec4 color;
	float t;
};

//============
//node information utility
//Utility
nodeTexel intTexture3D(usampler3D tex, uint z, uint y, uint x, uint maxZ, uint maxY, uint maxX)
{
	//maxZ (1:maxZ) number of element = maxZ
	// z = (0: max-1)
	float indexX = (float(x)+0.5/float(maxX))/float(maxX);
	float indexY = (float(y)+0.5/float(maxY))/float(maxY);
	float indexZ = (float(z)+0.5/float(maxZ))/float(maxZ);//TODO: convert to a vec operation

	uvec4 uSample = texture3D(tex, vec3(indexX, indexY, indexZ));

	nodeTexel node;
	node.first = uSample.r;
	node.second = uSample.a;
	return node;
}

uvec3 getNodeClusterPointer( nodeTexel node)
{
	uvec3 cPointer=uvec3(0u,0u,0u);
	cPointer[0] = node.first & uint(0x3FF) ;
	cPointer[1] = (node.first>>10u) & uint(0x3FF);
	cPointer[2] = (node.first>>20u) & uint(0x3FF);
	return cPointer;
}

nodeTexel getChildrenNode(uvec3 clusterPointer, uint k, uint j, uint i)
{
	nodeTexel node = intTexture3D(nodePool3D, clusterPointer[2]*2u+k, clusterPointer[1]*2u+j, clusterPointer[0]*2u+i,
													nodePoolDim,nodePoolDim,nodePoolDim);

	return node;
}

uvec3 getBlockPointer( nodeTexel node) 
{
	uvec3 bPointer=uvec3(0u,0u,0u);
	bPointer[0] = node.second & uint(0x3FF) ;
	bPointer[1] = uint( (node.second >> 10u) & uint(0x3FF));
	bPointer[2] = uint( (node.second >> 20u) & uint(0x3FF));
	return bPointer;
}

bool hasBlock( nodeTexel node)
{
	bool block = ((node.first>>30u &  uint(0x1)) >  uint(0))? true:false;
	return block;
}

bool isSubDividable( nodeTexel node)
{
	bool dividable =  ((node.first>>31u &  uint(0x1)) >  uint(0) )? true:false;
	return dividable;
}

//=======================
//Ray casting utility TODO:convert to ray representation
range intersectBox(vec3 entryPoint , vec3 rayDir, vec3 boxmin, vec3 boxmax)//, float *tnear, float *tfar)
{
    // compute intersection of ray with all six bbox planes
    vec3 invR = vec3(1.0) / rayDir;
    vec3 tbot = invR * (boxmin - entryPoint);
    vec3 ttop = invR * (boxmax - entryPoint);

    // re-order intersections to find smallest and largest on each axis
    vec3 tmin = min(ttop, tbot);
    vec3 tmax = max(ttop, tbot);

    // find the largest tmin and the smallest tmax
    float largest_tmin = max(max(tmin.x, tmin.y), tmin.z);
    float smallest_tmax = min(min(tmax.x, tmax.y), tmax.z);

	range tMinMax;
	tMinMax.tNear = largest_tmin;
	tMinMax.tFar = smallest_tmax;
	return tMinMax;
}

BlockBox findCurrentBlockAndBoundingBox_MultiRes(vec3 pos, float vt, float t, vec3 O) 
//vt is the distance between eye to the image plane, and t is the distance between eye to the current Position
{
	vec3 boxmin = vec3(0,0,0);
	vec3 boxmax;// = vec3(1,1,1);
	float boxsize = 1.0;//vec3(1,1,1);

	uvec3 nodeClusterPointer=uvec3(0u,0u,0u);
	uvec3 childrenIndex;

	nodeTexel currentNode;
	currentNode.first = 1u<<31u; currentNode.second = 0u;

	
	float scaleFactor = 1.0;
	
	bool SubD = true;
	bool isResOK = false;
	float voxelSize = 1.0/float(blockSize);
	BlockBox renderBlock;
	renderBlock.sampleRate = stepsize;
	renderBlock.LODSamplingFactor = 1.0;

	//while(SubD)
	while(SubD && !isResOK) //@@can't write isSubDividable(currentNode);@@
	{
		//decending
		childrenIndex.x = uint(int(2.0*pos.x));
		childrenIndex.y = uint(int(2.0*pos.y));
		childrenIndex.z = uint(int(2.0*pos.z));

		//check if current resolution is enough (the projected voxel size is smaller than a pixel in image space)
		if(1.0/length(boxmin - O)*voxelSize</*thresholdScreenPixelSize*/var)
			isResOK = true;
		voxelSize *= 0.5;

		boxsize *= 0.5;
		boxmin += vec3(childrenIndex)*boxsize;

		pos = abs(fract(2.0*pos));
		
		currentNode=getChildrenNode( getNodeClusterPointer(currentNode), childrenIndex.z, childrenIndex.y, childrenIndex.x);
		SubD = isSubDividable( currentNode );

		renderBlock.sampleRate *= 0.5;
		renderBlock.LODSamplingFactor *= 2;
	}

	uvec3 blockPointer=getBlockPointer(currentNode);
	boxmax = boxmin+vec3(boxsize);

	renderBlock.singleColor = float(currentNode.second)/voxelMax;  //due to the use of integer texture, the value is not in [0-1]
	renderBlock.hasBlock=hasBlock(currentNode);
	renderBlock.boxMin = boxmin;
	renderBlock.boxMax = boxmax;
	renderBlock.blockPointer = blockPointer;
	
	return renderBlock;
}

sample accumulateSingleColorBlock_showOctree(sample s,BlockBox renderBlock,range r)
{
	//accumlate
	vec4 value, texValue;
	vec3 pos;	
	texValue = texture1D(tranf, renderBlock.singleColor);
	texValue.a=0.1;
	//analytical integrate
	float d=(r.tFar - s.t)/stepsize;
	float t = (1.0-texValue.a);
	float T = exp(-t*d);
	s.color.rgb += T*texValue.rgb;
	s.color.a += T;

	return s;
}

 sample accumulateBlockRayVisualization(vec3 rayO, vec3 rayDir, sample s, BlockBox renderBlock, range r)
 {
	 s.t = r.tFar;
	 s.color.rgb =   (rayO+rayDir*r.tNear);
	 s.color.a = 1.0;
	 //s.color.rgb += (1.0-s.color.a)*0.03*(rayO+rayDir*r.tNear);
	 //s.color.a += (1.0-s.color.a)*0.03;
	 return s;
 }

 sample accumulateBlock_R(vec3 rayO, vec3 rayDir, sample s,BlockBox renderBlock, range r)
{
	//warping the coordinate
	uvec3 blockPointer = renderBlock.blockPointer;

	float invboxsize =1.0/( renderBlock.boxMax.r - renderBlock.boxMin.r);
	vec3 offset = (invBlockPoolDim)*vec3(blockPointer);
	float shrinkFactor = (float(blockSize-1u))/(float(blockSize));
	vec3 shrinkOffset = vec3(0.5/float(blockSize))*invBlockPoolDim;

	//accumlate
	float tEnd = r.tFar;
	vec4 value, texValue;
	vec3 pos;	

	int i=0;
	while(s.color.a < 0.99 && s.t < tEnd)
	{
		//safeguard
		//i++;
		//if(i>400)
		//{
		//	s.color = vec4(1,0,0,1);
		//	return s;
		//}

		pos = (rayO+s.t*rayDir);
		//return vec4(pos,1);
		
		pos = (pos-renderBlock.boxMin)*invboxsize;
		
		//return offset + (pos*invBlockPoolDim*shrinkFactor + shrinkOffset);
		
		//value = texture3D(blockPool, offset + (pos*invBlockPoolDim*shrinkFactor ));
		value = texture3D(blockPool, offset + pos*shrinkFactor*invBlockPoolDim + shrinkOffset);
		texValue = texture1D(tranf, value.a);
		//opacity correction   !!!!!!!!!!!!
		texValue.a = 1.0 - pow(1.0-texValue.a, renderBlock.sampleRate*1000);

		s.color.rgb += (1.0-s.color.a)*texValue.a*texValue.rgb;
		s.color.a += (1.0-s.color.a)*texValue.a;
		s.t += renderBlock.sampleRate;

	}

	//visualize block Pointer
	//color += vec4(offset,1)*0.2;

	return s;
}
  //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 sample accumulateBlockCount(sample s, BlockBox renderBlock, range r)
 {
	 s.color.rgb += (1.0-s.color.a)*0.05*vec3(1);//texValue.rgb;
	 s.color.a += (1.0-s.color.a)*0.05;
	 s.t += (floor((r.tFar - s.t)/renderBlock.sampleRate)+1)*renderBlock.sampleRate;
	 return s;
 }

sample accumulateBlockCountColor(vec3 rayO, vec3 rayDir, sample s, BlockBox renderBlock, range r)
 {
	//warping the coordinate
	uvec3 blockPointer = renderBlock.blockPointer;

	float invboxsize =1.0/( renderBlock.boxMax.r - renderBlock.boxMin.r);
	vec3 offset = (invBlockPoolDim)*vec3(blockPointer);
	float shrinkFactor = (float(blockSize-1u))/(float(blockSize));
	vec3 shrinkOffset = vec3(0.5/float(blockSize))*invBlockPoolDim;
	
	vec3 pos = (rayO+s.t*rayDir);
	//pos = (pos-renderBlock.boxMin)*invboxsize;
	//vec3 blockC = offset + pos*shrinkFactor*invBlockPoolDim + shrinkOffset;

	vec3 blockC = pos;
	blockC = vec3(1);
	float A = 0.02;
	s.color.rgb += (1.0-s.color.a)*A*blockC;//texValue.rgb;
	s.color.a += (1.0-s.color.a)*A;

	//this is why s.t += (floor((r.tFar - s.t)/renderBlock.sampleRate)+1)*renderBlock.sampleRate won't work
	//r.tFar can be smaller than s.t casue not ray increment 
	//if(r.tFar-s.t>0.0)
	//	s.t += (floor((r.tFar - s.t)/renderBlock.sampleRate)+1)*renderBlock.sampleRate; //less artifact
	//else
	//	s.t += 0.01*renderBlock.sampleRate;
	s.t += (floor(abs(r.tFar - s.t)/renderBlock.sampleRate)+1)*renderBlock.sampleRate;//renderBlock.sampleRate;

	return s;
 }
 
sample accumulateSingleColorBlock_Analytical(sample s,BlockBox renderBlock,range r)
{
	//accumlate
	vec4 value, texValue;
	vec3 pos;	
	texValue = texture1D(tranf, renderBlock.singleColor);

	//analytical integrate
	//TODO the analytical parameter is different 
	float F = exp(-texValue.a*(r.tFar - s.t)*240);
	s.color.rgb += (1.0-s.color.a)*(1.0 - F)*texValue.rgb;
	s.color.a += (1.0-s.color.a)*(1.0 - F);

	//r.tFar may smaller than s.t
	s.t += (floor(abs(r.tFar - s.t)/renderBlock.sampleRate)+1)*renderBlock.sampleRate;//renderBlock.sampleRate;
	return s;
}

sample accumulateBlock_AR(vec3 rayO, vec3 rayDir, sample s,BlockBox renderBlock, range r)
{
	//warping the coordinate
	uvec3 blockPointer = renderBlock.blockPointer;

	float invboxsize =1.0/( renderBlock.boxMax.r - renderBlock.boxMin.r);
	vec3 offset = (invBlockPoolDim)*vec3(blockPointer);
	float shrinkFactor = (float(blockSize-1u))/(float(blockSize));
	vec3 shrinkOffset = vec3(0.5/float(blockSize))*invBlockPoolDim;

	//accumlate
	float tEnd = r.tFar;
	vec4 value, texValue;
	vec3 pos;	

	int i=0;
	if(tEnd < s.t)
		tEnd += s.t-tEnd+0.01*renderBlock.sampleRate;
	while(s.color.a < 0.99 && s.t < tEnd)
	{
		//safeguard

		pos = (rayO+s.t*rayDir);
		//return vec4(pos,1);
		
		pos = (pos-renderBlock.boxMin)*invboxsize;
		
		value = texture3D(blockPool, offset + pos*shrinkFactor*invBlockPoolDim + shrinkOffset);
		texValue = texture1D(tranf, value.a);

		float F = exp(-texValue.a*renderBlock.sampleRate*renderBlock.LODSamplingFactor*30);
		s.color.rgb += (1-s.color.a)*(1.0 - F)*texValue.rgb;
		s.color.a += (1.0-s.color.a)*(1.0-F);
		s.t += renderBlock.sampleRate;

	}
	s.t += 0.01*renderBlock.sampleRate;
	//visualize block Pointer
	//color += vec4(offset,1)*0.2;

	return s;
}
 //===================================================================
vec4 visualizeOctree(vec3 start, vec3 end, vec3 O, vec3 screenNorm)
{
	vec4 color;
	float viewt = (1.0/dot(screenNorm,normalize(start-O)));
	float t = length(start-O);
	//viewt/t

	color = vec4(vec3(viewt/t),1);
///////////////////////////////////////////////////////
	vec3 rayDir = normalize(end - start);
	float tMax = length(end - start);
	
	int i=0;
	vec3 rayO = start;
	BlockBox renderBlock;
	range r ;
	vec3 pos;

	sample s;
	s.t = 15*0.0014;//stepsize;
	s.color = vec4(0,0,0,0);

	while(s.t<tMax-0.01 && s.color.a <0.99)
	{
		//ray decending, find the bounding box for current block
		pos = rayO + s.t*rayDir;
		renderBlock = findCurrentBlockAndBoundingBox_MultiRes(pos, viewt, t+s.t, O);
		//calculate ray box intersection

		r = intersectBox(rayO , rayDir, renderBlock.boxMin, renderBlock.boxMax);

		if(i==40)
		{
			//return vec4(1,0,0,1);
			break;
		}
		i++;

		//accumulate volume
		if(renderBlock.hasBlock)
		{
			//color = vec4(rayO+r.tNear*rayDir,1);
			//color = accumulateBlock(rayO, rayDir, color,renderBlock, r);
			//s = accumulateBlock_R(rayO, rayDir, s,renderBlock, r);
			s = accumulateSingleColorBlock_showOctree(s,renderBlock,r);
		}
		else
		{
			s = accumulateSingleColorBlock_showOctree(s,renderBlock,r);
		}

	}

	return s.color;
	//return color;
}
//////////////////////////////////////////////////////////////////
vec4 multiResOctree_V2(vec3 O, vec3 rayDir, float tStart, float tEnd)
{
	vec4 color;

	return color;
}

vec4 multiResOctree(vec3 start, vec3 end, vec3 O, vec3 screenNorm)
{
	vec4 color;
	float viewt = (1.0/dot(screenNorm,normalize(start-O)));
	float t = length(start-O);
	//viewt/t

	color = vec4(vec3(viewt/t),1);
	//return vec4(screenNorm,1);
	vec3 rayDir = normalize(end - start);
	float tMax = length(end - start);
	
	int i=0;
	vec3 rayO = start;
	BlockBox renderBlock;
	range r ;
	vec3 pos;

	sample s;
	s.t = 0.001*stepsize;//15*0.0014;//stepsize;
	s.color = vec4(0,0,0,0);

	while(s.t<tMax/*-0.01*/ && s.color.a <0.99)
	{
		//ray decending, find the bounding box for current block
		pos = rayO + s.t*rayDir;
		renderBlock = findCurrentBlockAndBoundingBox_MultiRes(pos, viewt, t+s.t, O);
		//calculate ray box intersection

		r = intersectBox(rayO , rayDir, renderBlock.boxMin, renderBlock.boxMax);

		//accumulate volume
		if(renderBlock.hasBlock)
		{
			//color = vec4(rayO+r.tNear*rayDir,1);
			//color = accumulateBlock(rayO, rayDir, color,renderBlock, r);
			//if(integralMethod>0.5)
			//	s = accumulateBlock_R(rayO, rayDir, s,renderBlock, r);				    
			//else
			//s = accumulateBlockCountColor(rayO, rayDir, s, renderBlock,r);
			s = accumulateBlock_AR(rayO, rayDir, s,renderBlock, r);	//accurate integration
			//s = accumulateBlockRayVisualization(rayO, rayDir, s, renderBlock, r);
			//return s.color;
		}
		else
		{
			//s = accumulateBlockCountColor(rayO, rayDir, s, renderBlock,r);
			s = accumulateSingleColorBlock_Analytical(s,renderBlock,r);
			//s = accumulateBlockRayVisualization(rayO, rayDir, s, renderBlock, r);
			//s.t = r.tFar;
			//s = accumulateSingleColorBlock_showOctree(s,renderBlock,r);
			//return s.color;
		}

	}

	return s.color;
	//return color;
}

void main()
{
	vec3 boxMin = vec3(-1.0, -1.0, -1.0);
    vec3 boxMax = vec3(1.0, 1.0, 1.0);

	//inside volume view
	mat3 rotateMat = mat3(vec3(invViewMatrix[0]),vec3(invViewMatrix[1]),vec3(invViewMatrix[2]));
	vec4 O;
	if(insideView>0.5)
	{
		vec3 rotTranslate = rotateMat*translateVec3;
		O = (vec4(-rotTranslate[0],-rotTranslate[1],4-rotTranslate[2],1));
	}
	else
	{
	//normal view
		O = invViewMatrix*(vec4(-translateVec3[0],-translateVec3[1],4-translateVec3[2],1));
	}
	vec3 eyeO = O.xyz/O.w;


	vec2 uv = (gl_TexCoord[0].xy/gl_TexCoord[0].w);
	uv[0] *= aspectRatio;
	vec3 eyeD = rotateMat*normalize(vec3(uv,-2.0));

	vec3 screenNorm = rotateMat*normalize(vec3(0,0,-2.0));
	
	range t = intersectBox(eyeO, eyeD, boxMin, boxMax);

	vec3 start = vec3(0);
	vec3 end = vec3(0);
	
	if(t.tNear<0)
		t.tNear = 0;
	
	if(t.tNear>t.tFar)
		gl_FragColor = vec4(0,0,0,1);
	else
	{
		//gl_FragColor = vec4(pos*0.5+0.5,t.tNear);
		start = (eyeO + eyeD*t.tNear)*0.5+0.5;
		end = (eyeO + eyeD*t.tFar)*0.5+0.5;
		vec3 correctedO = eyeO*0.5+0.5;

		ray r;
		r.o = eyeO; r.d = eyeD;

		if(isoValue > 0.5)
			gl_FragColor = visualizeOctree(start, end, correctedO,screenNorm);
		else
			gl_FragColor = multiResOctree(start, end, correctedO,screenNorm);
	}

}